#include <stdio.h>
#include "common.c"

unsigned int phred_for_delete = 15;
//first 3 input: CIGAR="3M2I2M2D2M", only M I D and I followed by M, number(even 0) between these char; last 2 output.
unsigned int sites_allele_quality(const unsigned char *cigar, const unsigned char *read_seq, const unsigned char *base_quality, unsigned char allele_seq[][ALLELE_LEN], unsigned int allele_quality[]) {
    unsigned int num_site = 0;//next site index.

    unsigned int prev_insert_num = 0;//!!indicate
    unsigned char prev_insert_seq[ALLELE_LEN - 2];//no match base and \0.
    unsigned int prev_insert_quality;
    unsigned long int num;
    unsigned char mdi;
    while(num = strtol(cigar, &cigar, 10), (mdi = *cigar ++) != '\0') {
        //warning, strtol 2nd arg expect "char **"
        if(mdi == 'M') {
            unsigned int i;
            for(i = 0; i < num; i ++, num_site ++) {
                if(i == 0 && prev_insert_num > 0) {
                    strncpy(allele_seq[num_site], prev_insert_seq, prev_insert_num);
                    allele_seq[num_site][prev_insert_num] = *read_seq ++;
                    allele_seq[num_site][prev_insert_num + 1] = '\0';
                    allele_quality[num_site] = prev_insert_quality + (*base_quality ++ - 33);
                    prev_insert_num = 0;//tell next match base !!
                } else {
                    allele_seq[num_site][0] = *read_seq ++;
                    allele_seq[num_site][1] = '\0';
                    allele_quality[num_site] = (*base_quality ++ - 33);
                }
            }
            prev_insert_num = 0;//no need.
        } else if(mdi == 'I') {
            prev_insert_num = num;//!!
            strncpy(prev_insert_seq, read_seq, num);
            read_seq += num;
            prev_insert_quality = 0;//!!
            for(; num > 0; num --)
                prev_insert_quality += (*base_quality ++ - 33);
        } else if(mdi == 'D') {
            for(; num > 0; num --, num_site ++) {
                allele_seq[num_site][0] = '*';
                allele_seq[num_site][1] = '\0';
                allele_quality[num_site] = phred_for_delete;
            }
            prev_insert_num = 0;
        } else {
            fprintf(stdout, "CIGAR: %c\n", mdi);
            break;//no
        }
    }
    return num_site;
}//OK!! "2M2I2M2D2M", "ATCCGTGG", "########"

void add_allele(struct sample_site *site, const unsigned char *seq, const unsigned int quality) {
    unsigned int i;
    for(i = 0; i < site->num_allele && strcmp(site->allele_seq[i], seq) != 0; i ++) ;
    if(i < site->num_allele) {
        site->allele_reads[i] += 1;
        site->allele_quality[i] += quality;
    } else {
        //now i = num_allele.
        strcpy(site->allele_seq[i], seq);
        site->allele_reads[i] = 1;//initial
        site->allele_quality[i] = quality;
        site->num_allele ++;
     }
}
void add_allele_to_sites(unsigned int num_site, const unsigned char allele_seq[][ALLELE_LEN], const unsigned int allele_quality[], struct sample_site *from) {
    unsigned int i;
    struct sample_site *prev;//first from_site not null !
    for(i = 0; i < num_site; i ++) {
        if(from == NULL) {
            from = malloc( sizeof(struct sample_site) );
            from->chr = prev->chr;
            from->pos = prev->pos + 1;
            from->num_allele = 0;
            from->next = NULL;
            prev->next = from;
        }
        add_allele(from, allele_seq[i], allele_quality[i]);
        prev = from;
        from = from->next;
    }
}

void output_site_divide(unsigned int chr, unsigned int pos, unsigned int num, const unsigned char seq[][ALLELE_LEN], const unsigned int reads[], const unsigned int quality[]) {
    fprintf(stdout, "%u\t%u\t", chr, pos);
    unsigned int i;
    for(i = 0; i < num; i ++)
        fprintf(stdout, "%s:%u:%u/", seq[i], reads[i], quality[i] / strlen(seq[i]) );
    fprintf(stdout, "\n");
}


int main(int argc, char *argv[]) {
    unsigned int min_map_quality = atoi( argv[1] );
    phred_for_delete = atoi( argv[2] );

    struct sample_site *from = NULL;
    #define LINE_LEN 200000
    unsigned char line[LINE_LEN];

    while(fgets(line, LINE_LEN, stdin)) {
        unsigned int flags;
        unsigned char chr_str[6];
        unsigned int pos, map_quality;
        #define CIGAR_LEN 10000
        unsigned char cigar[CIGAR_LEN];
        #define READ_LEN 2000
        unsigned char read_seq[READ_LEN];
        unsigned char base_quality[READ_LEN];
        sscanf(line, "%*s\t%d\t%s\t%d\t%d\t%s\t%*s\t%*d\t%*d\t%s\t%s\t", &flags, chr_str, &pos, &map_quality, cigar, read_seq, base_quality);

        if((flags & 4) || (map_quality < min_map_quality))
            continue;
        unsigned int chr = get_chr_num(chr_str);//!!

        for(; (from != NULL) && compare_pos(from->chr, from->pos, chr, pos) < 0; ) {
            output_site_divide(from->chr, from->pos, from->num_allele, from->allele_seq, from->allele_reads, from->allele_quality);
            struct sample_site *save = from->next;
            free(from);
            from = save;
        }
        if(from == NULL) {
            from = malloc( sizeof(struct sample_site) );
            from->chr = chr;
            from->pos = pos;
            from->num_allele = 0;
            from->next = NULL;
        }//now from is chr & pos...

        #define READ_SITES_NUM 500
        unsigned char allele_seq[READ_SITES_NUM][ALLELE_LEN];
        unsigned int allele_quality[READ_SITES_NUM];
        unsigned int num_site = sites_allele_quality(cigar, read_seq, base_quality, allele_seq, allele_quality);
        add_allele_to_sites(num_site, allele_seq, allele_quality, from);
    }
    
    for(; from != NULL; ) {
        output_site_divide(from->chr, from->pos, from->num_allele, from->allele_seq, from->allele_reads, from->allele_quality);
        struct sample_site *save = from->next;
        free(from);
        from = save;
    }
}
